/* Histogram for N-D arrays 

   [X, N] = NDHIST(Y, WIDTH) sorts the data in Y into equally-
   spaced bins of size WIDTH, returning the bin values in X and
   the bin counts in N.  Y can have arbitrary shape & size.  X
   and N are returned as DOUBLE, regardless of Y.
   
   [X, N] = NDHIST(Y, WIDTH, RANGE) generates bins over the
   range specified in the two element vector RANGE.

/* Ignores imaginary parts of complex inputs */

/* By:   S.C. Molitor (smolitor@bme.jhu.edu)
   Date: February 26, 1999 */


// inclusions

#include "mex.h"
#include "uintutil.h"


// MEXFUNCTION - gateway routine

void mexFunction(int			nlhs,
				 mxArray		*plhs[],
				 int			nrhs,
				 const mxArray	*prhs[]) {

	// declare variables

	mxArray		*Y, *W, *R;		// pointers to input data, width & range
	mxArray		*X, *N;			// pointers to output bin values & counts
	mxClassID	classID;		// input Y class ID
	double		*prY;			// pointer to input data
	double		*prX, *prN;		// pointers to bin values & counts
	double		bwidth;			// bin value increment
	double		bmin, bmax;		// min & max bin values
	double		tval;			// temporary value
	int			lengthY;		// length of input data
	int			lengthN;		// number of bins
	int			indexN;			// index of N to add count
	int			i;				// loop counter	

	// call VALIDATEINPUTS to validate calling syntax

	ValidateInputs(nlhs, nrhs, prhs, 2, 3, 2);

	// determine length of input array Y
	// if empty, return as empty outputs

	Y = (mxArray *) prhs[0];
	lengthY = mxGetNumberOfElements(Y);
	if (lengthY <= 0) {
		plhs[0] = mxDuplicateArray(Y);
		plhs[1] = mxDuplicateArray(Y);
		return;
	}

	// obtain pointer to input Y if DOUBLE
	
	classID = mxGetClassID(Y);
	if (classID == mxDOUBLE_CLASS) {
		prY = (double *) mxGetData(Y);
	}
	
	// otherwise create a DOUBLE input
	// convert each point individually
	
	else {
		prY = (double *) mxCalloc(lengthY, sizeof(double));
		for (i = 0; i < lengthY; i++) {
			prY[i] = GetDouble(Y, classID, i);
		}
	}

	// obtain bin width

	W = (mxArray *) prhs[1];
	if (mxGetNumberOfElements(W) == 1) {
		bwidth = GetDouble(W, mxGetClassID(W), 0);
		if (bwidth <= 0.0) {
			mexErrMsgTxt("WIDTH must be positive");
		}
	}
	else {
		mexErrMsgTxt("WIDTH must be scalar");
	}

	// obtain bin range if provided

	if (nrhs > 2) {
		R = (mxArray *) prhs[2];
		if (mxGetNumberOfElements(R) == 2) {
			bmin = GetDouble(R, mxGetClassID(R), 0);
			bmax = GetDouble(R, mxGetClassID(R), 1);
			if (bmin > bmax) {
				tval = bmin;
				bmin = bmax;
				bmax = tval;
			}
		}
		else {
			mexErrMsgTxt("RANGE must be two element vector");
		}
	}

	// otherwise find range the old-fashioned way

	else {
		bmin = prY[0];
		bmax = prY[0];
		for (i = 1; i < lengthY; i++) {
			tval = (double) prY[i];
			if (tval < bmin) {
				bmin = tval;
			}
			else if (tval > bmax) {
				bmax = tval;
			}
		}
	}

	// determine number of bins
	// create output arrays

	lengthN = (int) ((bmax - bmin) / bwidth) + 1;
	X = mxCreateDoubleMatrix(1, lengthN, mxREAL);
	plhs[0] = X;
	prX = (double *) mxGetData(X);
	N = mxCreateDoubleMatrix(1, lengthN, mxREAL);
	plhs[1] = N;
	prN = (double *) mxGetData(N);

	// create bin values

	tval = bmin;
	for (i = 0; i < lengthN; i++) {
		prX[i] = tval;
		tval += bwidth;
	}

	// sort data in bins to obtain bin counts
	// data is offset by MIN & normalized by WIDTH to create index!
	// data outside the range of bins is ignored!

	for (i = 0; i < lengthY; i++) {
		indexN = (int) (((double) prY[i] - bmin) / bwidth);
		if ((indexN >= 0) && (indexN < lengthN)) {
			prN[indexN] += 1.0;
		}
	}

	// free pointer to Y if one was created

	if (classID != mxDOUBLE_CLASS) {
		mxFree((void *) prY);
	}
}
